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1. Introduction 

This paper considers the problem of solving the Helmholtz equation 

- /\il){x,y) = Eip{x,y) (1) 

over a two-dimensional domain, of arbitrary shape, assuming Dirichlet boundary 
conditions over the border, OB. Physically, this equation describes the classical vibration 
of a homogenoeous membrane or the behaviour of a particle confined in a region with 
infinite walls in quantum mechanics. Unfortunately exact solutions to this equation are 
available only in few cases, such as for a rectangular or a circular membrane, where 
they can be expressed in terms of trigonometric and Bessel functions respectively pQ. 
In the majority of cases, in fact, only numerical approaches can be used: some of these 
approaches are discussed for example in a beautiful paper by Kuttler and Sigillito, [2]. 
The purpose of the present paper is to introduce a different approach to the numerical 
solution of the Helmholtz equation (both homogenous and inhomogeneous) and illustrate 
its strength and flexibility by applying it to a large number of examples. 

The paper is organized as follows: in Section |2] I describe the method and discuss 
its application to the classical problem of a L-shaped membrane; in Section |3] I consider 
an homogenous membrane, with the shape of Africa and calculate few states; in Section 
|4] I consider two inequivalent membranes, which are known to be isospectral, obtaining 
a numerical indication of isospectrality; in Section [5] I study an example of irregular 
drum; in Section |6] the method is applied to study the emergence of bound states in 
a configuration of wires of neglegible transverse dimension, in presence of crossings; 
in Section [7] I show that even more precise results can be achieved by combining the 
collocation method with a conformal mapping of the boundary. Finally, in Section |8] I 
draw my conclusions. 

2. The method 

The method that I propose in this paper uses a particular set of functions, the Little 
Sine functions (LSF) of [13l HI], to obtain a discretization of a finite region of the 
two-dimensional plane. These functions have been used with success in the numerical 
solution of the Schrodinger equation in one dimension, both for problems restricted 
to finite intervals and for problems on the real line. In particular it has been proved 
that exponential convergence to the exact solution can be reached when variational 
considerations are made (see [131 E]). 

Although Ref.[T3] contains a detailed discussion of the LSF, I will briefly review 
here the main properties, which will be useful in the paper. Throughout the paper I 
will follow the notation of [13]. 

A Little Sine Function is obtained as an approximate representation of the Dirac 
delta function in terms of the wave functions of a particle in a box (being 2L the size 
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of the box). Straightforward algebra leads to the expression 

2N y sinx_(a;) cosx+(a^) J 

where x±{.^) = 2M^(^ ^ ^^)- The index k takes the integer values between —N/2 + 1 
and N/2 — 1 (A^ being an even integer). The LSF corresponding to a specific value of 
k is peaked at Xk = 2Lk/N = kh, h being the grid spacing and 2L the total extension 
of the interval where the function is defined. By direct inspection of eq. ^ it is found 
that Sk{h,N,Xj) = Skj, showing that the LSF takes its maximum value at the k^'^ grid 
point and vanishes on the remaining points of the grid. 

It can be easily proved that the different LSF corresponding to the same set are 
orthogonal [13]: 

Sk{h,N,x)sj{h,N,x)dx = h 6kj (3) 



-L 

and that a function defined on a; G {—L, L) may be approximated as 

N/2-l 

f{^)^ E fi^k) Sk{h,N,x) . (4) 

k=-N/2+l 

This formula can be applied to obtain a representation of the derivative of a LSF 
in terms of the set of LSF as: 

dsk{h,N,x) dsk{h, N,x) 



d'^Sk{h, N,x) ^ d'^Sk{h, N,x) 



dx"^ 

j 



s,{h,N,x)^Y.''kj sj{h,N,x) 



(5) 

where the expressions for the coefficients c^i^j can be found in [13]. Although eqs.(4) is 
approximate and the LSF strictly speaking do not form a basis, the error made with 
this approximation decreases with and tends to zero as tends to infinity, as shown 
in [13]. For this reason, the effect of this approximation is essentially to replace the 
continuum of a interval of size 2L on the real line with a discrete set of — 1 points, 
Xk, uniformly spaced on this interval. 

Clearly these relations are easily generalized to functions of two or more variables. 
Since the focus of this paper is on two dimensional membranes, I will briefly discuss how 
the LSF are used to discretize a region of the plane; the extension to higher dimensional 
spaces is straightforward. A function of two variables can be approximated in terms of 
{N^ — 1) X [Ny — 1) functions, corresponding to the direct product of the A^^,' ~ 1 and 
Ny — 1 LSF in the x and y axis: each term in this set corresponds to a speciflc point 
on a rectangular grid with spacings hx and hy (in this paper I use a square grid with 
Nx = Ny = N and = Ly = L). 
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Since (/c, k') identifies a unique point on the grid, I can select this point using a 
single index 



f N 

K = k' + — + {N -l)\k+ 1 



(6) 



which can take the values 1 < K < [N 

K 

k = l-N/2 + 



lY. I can also invert this relation and write 



k' 



_N -l + e 
K - N/2 - {N - 1) 



K 



(7) 
(8) 



0. 



_N-l + e_ 

where [a] is the integer part of a real number a and e 

As a natural extension of the results presented in [131 113] I can consider the 
Schrodinger equation in two dimensions 

Hipnix, y) = [-A + V{x, y)] ipnix, y) = Enipnix, y) (9) 

using the convention of assuming a particle of mass m = 1/2 and setting h = 1. The 
Helmholtz equation, which describes the vibration of a membrane, is a special case of 
(|9|, corresponding to having V{x^y) = inside the region B where the membrane lies 
and V{x,y) = oo on the border dB and outside the membrane. 

The discretization of eq. ^ proceeds in a simple way using the properties discussed 
m eqs. Q and (§: 



kk',jj' — ~ Ckj^k'j'+SkjC)^iji + SkjSk'j'V{Xk,yk') (10) 

where {k,j,k',j') = —N/2 + l,...,N/2 — 1. Notice that the potential part of the 
Hamiltonian is obtained by simply "collocating" the potential V{x,y) on the grid, an 
operation with a limited computational price. The result shown in (10) corresponds to 
the matrix element of the Hamiltonian operator H between two grid points, {k, k') and 
(j, j'), which can be selected using two integer values K and J, as shown in ([6]). 

Following this procedure the solution of the Schrodinger (Helmholtz) equation 
on the uniform grid generated by the LSF corresponds to the diagonalization of a 
(A^ — 1)^ X (A^ — 1)^ square matrix, whose elements are given by eq. (10). 

I will now use a specific problem, the vibration of a L-shaped membrane, represented 
in FigjT| to illustrate the method, and discuss different implementations of the method 
itself. This problem has been widely used in the past to test the performance of the 
different numerical methods (see for example refs. [SI HI |5|, [21 [71 13 El [TOl [11]) and it is 
therefore a useful tool to assess the strength of the present approach. Because of the 
reentrant corner, corresponding to the angle = 37r/2 located at (0,0), the derivatives 
of ip{x,y) in the radial direction are unbounded (see [3]). 

Reid and Walsh in [3] obtained a numerical approximation for the two lowest modes 
of this membrane using finite differences and a confomal map which eliminates the 
reentrant corner (see fig. 5 of [3]); a more precise result was later obtained by Fox, Henrici 
and Moler who used the Method of Particular Solutions (MPS) in [1] exploiting the 
symmetries of the problem (the reader may find a detailed discussion of the symmetries 



.(2) 
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Figure 1. L shaped membrane. The dots are the collocation points corresponding to 
N = 14. 



for this problem in [2j): the first eight digits of the lowest eigenvalue reported by the 
authors are correct. Mason has obtained numerical estimates for the first few modes of 
the L-shaped membrane in terms of a two dimensional Chebyshev series [S]. Milsted 
and Hutchinson |6] have obtained finite element solutions to this problem. Sideridis in 
[7] used a conformal mapping of the L-shaped region onto a square and then solved the 
resulting equation on a uniform rectangular mesh, obtaining the first four digits of the 
lowest mode. Schiff, ref. [8], has calculated the first 15 lowest modes of this membrane 
using finite elements, with a refined grid covering the region surrounding the reentrant 
corner. 

More recently Platte and Driscoll have solved the boundary value problem on the 
L-shape membrane using radial basis functions [9]. Finally Betcke and Trefethen have 
revisited the MPS in [10]; in that paper they have observed that the MPS reaches 
a minimal error for a certain value of (the number of collocation points on each 
of the sides non adjacent to the corner where the expansion is performed) but then 
it starts to grow as increases. The modified version of the method discussed in 
[To] , which samples the Fourier-Bessel functions also in the interior points, corrects this 
problem and provides a convergent behaviour for the error. In this way Betcke and 
Trefethen were able to obtain the first 14 digits of the lowest eigenvalue of the L-shaped 
membrane, Ei k, 9.6397238440219. I will use this precise result to test the accuracy of 
our method. Ref.[TT] contains precise estimates for some higher excited states of the 
L-shaped membrane. 

I will now apply the LSF to the numerical solution of this problem: looking at 
Fig. [T] I consider the grid points which are internal to the membrane and which do not 
fall on the border. For a fixed there is a total of 3/4A^^ — 2A^ + 1 points; the grid 
represented in the figure corresponds to = 14 and therefore to a total of 120 internal 
points. In this case the collocation of the Hamiltonian on the uniform grid generated by 
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Table 1. First 108 eigenvalues of the L-shaped membrane calculated with a grid with 
N = 60. 
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the LSF leads to a 120 x 120 matrix, which can then be diagonahzed. The eigenvalues 
of this matrix provide the lowest 120 modes of the membrane, while the eigenvectors 
provide the lowest 120 wave functions. Alternatively I can pick all the points of the grid 
internal to the membrane, including those falling on the border: in such case a total of 
3/4iV^ — points is found, corresponding to a total of 133 points in the case of the 
Figure. 

Table [2] contains the first 108 eigenvalues of the L-shaped membrane calculated 
using a grid with = 60 and selecting the grid points according to the prescriptions 
just explained. I have used the notation E^"* for the energy of the n*^ state when the 
collocation points on the border are either rejected {e'^'') or kept (E'i"^). The notation 
(±) is used since the two sets approach the exact results either from above (+) or from 
below (— ), as one can see comparing these numbers with the precise results contained 
in [Tot [TT] . The reader will certainly notice that the results of Table [2] contain rather 
large errors: in the case of the fundamental state, for example, one has an error of about 
1% from E'i^^ and a much larger error of almost 5% for En ^ 

The left panel of Fig|2]shows the eigenvalues E^t^ (solid line) and i^n"-* (dashed line) 
for the L-shaped membrane corresponding to a grid with A^ = 60. The reader may notice 
that the higher end of the spectrum displays a curvature, contrary to the behaviour 
predicted by Weyl's law, i.e. (A^) oc E for large energies. It is easy to show that such 
effect is artificial: consider for example the case of a particle confined in a unit square, 
whose energies are given by En^ „ = (n^+?7,y)7r^. The diagonalization of the Hamiltonian 



\2 



(10) for this problem would provide the energies corresponding to the (A^ — 1) states 
obtained taking the first A^ — 1 values of Ux and riy. This means that for energies 
higher than E^q = [N"^ — 2N + 2] vr^ the method will provide only the eigenvalues 
contained inside a square of side A^ — 1 (in the (n^, ny) plane), up to a maximal energy 
Em AX = 2 [A^^ — 2A^ + 2] tt^. For this reason, the states above E^ are incomplete and 
should not be taken into account for inferring the asymptotic behavior of (A^). The right 
panel of Fig, 2 displays the asymmetry defined as An = 2(£'i"^^ — Eii~^)/{E^'^ + i?i~^) 
for the same grid: this quantity provides an upper estimate for the error. 

Fig. [3] displays the ground state energy of the L-shaped membrane as a function 
of the number of grid points and compares it with the precise result of [lOj: as already 
pointed out the two sets approach the exact value from above and below. 

Much more precise results can be obtained by performing an extrapolation of the 
results corresponding to finite grids: this is a common procedure used in the literature 
(see for example [2]). I have considered four different extrapolation sets using the 
numerical results obtained working with grids with A^ ranging from A^ = 10 to A^ = 60 
(only even values). Calling h = 2L/N the grid spacing the sets are: 

N 

f,{h) = Y,Cnh^ (11) 

n=0 

f2{h) = ^"=°;f ^ (12) 

1 + Eni? Cnh- 
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19.7392088020095 


o 

6 


(+) 


19.7392088019879 


19.7392088021704 


19.7392087962239 


T n Tononoonoi ^or* 

19.7392088021785 


A 

4 


(-) 


29.5178267971821 


29.5214811097206 




on r'O"! /ioi"i"ino/io 

29.5214811103487 


4 


(+) 


29.5214810813053 


29.5214811126514 


on r o "1 /I '7n a r ono 1 

29.5214794563921 


on roi /I o 1 "1 "1 /I "1 rn/^* 

29.5214811141506 


5 


(-) 


31.9159767579531 


31.9125745966885 




01 ni o/^op'nr'oonor'if: 

31.9126359533035 


5 


(+) 


31.9123209946513 


31.9126005580344 


31 .9126386707453 


OI ni o/^ornr^i o/^o* 

31.9126359571263 





(-) 


41.474267306813 


/IT A ''7 A /l^/innOOOT O 

41.4744740922213 


/II A 'in Tni /1/100000 

41.4761914432832 


A ^ /iT/irnnni /i o ^Tn* 

41.4745099148779 


6 


(+) 


41 .4742769974:4:02 


41.4744780007070 


41.4741677038785 


A -^ A 1 A r\r\c>f\r\ A A c> 

41.4745098904487 


on 




ini 77fic;Ri fivciQi /I* 

iUi.( (D00i0(0oi4 


ini finc;QQQQfioo7c; 
iUi.DUOoooooyy ( 




yy. ( ( iOZZ4o0iUoO 


20 


(+) 


101.604853531780 


101.605223692426 


101.673183488214 


101.605294080845* 


50 






246.740564791939 




246.602432808866* 


50 


(+) 


250.784799377301 


250.785244396338 




250.785494606618* 


104 






410.08260648211 






104 


(+) 


493.480067984180* 


493.480206216096 




493.488405725447 



Table 2. Extrapolation of the nine eigenvalues of the L-shaped membrane using the 
four different sets. The first 6 states correspond to extrapolating the results for grids 
going from N = 10 to = 60, with 25 unknown coefficients; the last two states 
correspond to extrapolating the results for grids going from = 18 to = 60, and 
with 21 unknown coefficients. For a given state, the set with the asterisk corresponds 
to the minimal value taken by the least squares. The results which do not converge to 
the exact value have been omitted. 



N 



fs{h)=Co + J2cnh-/'^'/' (13) 



n=l 



C„ + ^„/3+2/3 

where N is an even integer which determines the number of coefficients used in the fits. 
The continuum hmit is reached taking /i — > oo, where only the coefficient Cq 



survives. The unknown coefficients in the expressions (11), (12), (13) and (14) are 
obtained using a Least Square approach: I show the results of this procedure in Table 
[2j In general, the last set provides the best results and indeed it reproduces the ffist 11 
digits of El correctly, using either the values of e[ ^ or those of e['^\ In the case of E^, 
for which the exact value is known {E^ = 2-7?^), I obtain the ffist 14 digits correct using 
E^'^'' and the ffist 11 digits correct using E^ \ 
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Figure 2. Left pancLEnergy of the ground state of the L shaped membrane as a 
function of the number of grid points N . The horizontal hne is the precise result of 
[in]. The set approaching the exact result from above (below) corresponds to 
{e[~'^). Right panel: The asymmetry An = 2(i;i+^ - eI~^) / {eI^'^ + E^n^) calculated 
with a grid with N — 60. 
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Figure 3. Energy of the ground state of the L shaped membrane as a function of 
the number of grid points N . The horizontal line is the precise result of [lOj. The set 
approaching the exact result from above (below) corresponds to e'\^^ {e\ ■*). 



In [15] Berry has devised an algorithm for obtaning successive approximations to 
the geometric properties Kj of a closed boundary B given the lowest eigenvalues 
En- The partition function $(t) = e~^"* obeys an asymptotic expansion for small 

values of t 

<^.(t)^-5^ir,t^•/^ (15) 

i=o 

where the coefficients Kj are related to the geometric properties of B. For example 
= A/Att and Ki = —7^/8-^/7?. Using this asymptotic expansion Berry has obtained 
accelerated expressions for the geometrical constants of B. In particular for the area of 
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B he has found the approximant (eq.(20) of [15] ) 

m! ^ — ^ 

where ^„ = \/Ent. 
Aiit) 



n=l 



Llit) 



10 



(16) 




Figure 4. Left Panel: The area approximant A2{t) obtained using the expression of 
Berry. The thin solid and dashed lines are obtained with the first 1000 eigenvalues 
corresponding to the sets E^^'' and i?i ■* respectively. The bold solid and dashed lines 
correspond to the sets obtained through an extrapolation from the original sets. Right 
Panel: The perimeter approximant L2{t) obtained with the improved expression of 
Berry. The same sets of eigenvalues have been considered. 

In the left panel of Fig|4] 1 show the area approximant ^2(^)5 obtained using the 
expression of Berry. The thin lines correspond to using the the sets E^'^ and En ■* 
(solid and dashed lines respectively); the thick lines correspond to using the eigenvalues 
obtained from the extrapolation of the sets e'^^'^ and En ^ (solid and dashed lines 
respectively). 1 call E^^ the eigenvalues obtained extrapolating the eigenvalues E^"^-^ 
the extrapolation is carried out using the results obtained with grids with N going from 
= 48 to = 60 and assuming En{N) ^ En + ^ ■ The approximants obtained with 
the extrapolated eigenvalues provide excellent approximations to the area and perimeter 
of the membrane, as seen in Fig|4j 

Fig. [5] shows the first two eigenfunctions of the L shaped membrane obtained 
with a grid corresponding to = 30. The solid lines appearing in the "forbidden 
region" correspond to the level ip{x,y) = 0: the effect observed in the figure is due to 
the approximation of working with a finite number of grid points. In fact, although a 
particular LSF vanishes on the points defining the grid, except on a particular point, 
where it reaches its maximum, it is non-zero elsewhere. This means that the numerical 
solution can take small values even in the region where the exact solution must vanish; 
however, the size of this effect decreases as the number of grid points is increased (taking 
into account that the computational load roughly increases as A^*^). In the Appendix we 
propose an alternative procedure which does not involve the diagonalization of larger 
matrices and which can be used to improved the results obtained with a given grid. 
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3. The Africa drum 

I will now examine the case of a membrane with an irregular shape. The application of 
the method proceeds exactly as in the case of the L-shaped membrane: once a grid is 
chosen, the points of the grid which are internal to the membrane are used to build a 
matrix representation of the Hamiltonian which, once diagonalized, provides the energies 
and wave functions of the problem. 

As a paradigm of this class of membranes I have studied the vibrations of a drum 
with the shape of Africa. Unlike in the previous example the border does not cross 
the grid points, a feature which affects the precision of the results. The plots in Fig|6] 
display the energies of the first two states of the Africa drum for grids with different 
N (the dots in the plots) and compare them with the best fit obtained assuming that 
E{N) = a + b/N, where a and b are constants independent of N. The irregularity of 
the border is reflected in the behavior of the eigenvalues which decay with N but at the 
same time oscillate. 

In Fig|7]l show the density plot of four different states of the Africa drum, obtained 
using a grid with N = 60. In Fig|8]I show the wave function of the ground state of the 
Africa drum, obtained using a grid with N = 60. 

4. Isospectral membranes 

In a classic paper dated 1966, [16j, Kac formulated an interesting question: whether 
it is possible to hear the shape of a drum, meaning if the spectrum of frequencies of 
a given drum is unique to that drum or drums with different shapes can have the 
same spectrum. The question was flnally answered in 1992, when Gordon, Webb and 
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Figure 6. Left: Energy of the fundamental mode of the Africa shaped membrane as a 
function of the number of grid points. The continous Hne is the fit i?i = a + h/N , with 
a = 20.1705. Right: Energy of the first excited mode of the Africa shaped membrane 
as a function of the number of grid points. The continous line is the fit Ei — a + b/N, 
with a ^ 32.2774. 
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Figure 7. Density plot for the fundamental state (upper left), first excited state 
(upper right), 200*'* excited state (lower left) and 300*'' excited state (lower right) of 
the Africa shaped membrane. In all plots the absolute value of the wave function is 
shown and a grid with = 60 is used. 
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Figure 8. Ground state of the Africa shaped membrane obtained using = 60. 

Wolpert found a first example of inequivalent drums having the same spectrum [T7]. 
An experiment made by Sridhar and Kudrolh reported in [18] used microwave cavities 
with the shape of the drums of |17j to verify the equality of the spectrum for the lowest 
54 states. More recently the same experiments have been carried out on isospectral 
cavities where the classical dynamics changes from pseudointegrable to chaotic [19]. 
Numerical calculations of the first few modes of the isospectral drums found in [T7] have 
been performed with different techniques: Wu, Sprung and Martorell [20] have used a 
mode matching method to calculate the first 25 states of these drums and compared the 
results with those obtained with finite difference; using a different approach DriscoU [21] 
has also calculated the first 25 states obtaining results which are accurate to 12 digits; 
Betcke and Trefethen [TD] have used their modified version of the method of particular 
solutions to obtain the first three eigenvalues of these drums, reporting results which 
are slightly more precise than those of DriscoU. 

I will now discuss the application of the present method to the calculation of the 
spectrum of these isospectral membranes: whereas in the case of the L-shaped membrane 
the border of the membrane was sampled by the grid, regardless of the grid size (keeping 
even), in the case of the isospectral membranes this happens only for grids where 
A^ = 6fc, with k integer. It is important to restrict the calculation to this class of grids to 
avoid the oscillations observed in the case of the Africa membrane. I have thus applied 
the method with grids ranging from A^ = 6 to A^ = 12C[§j 

§ The numerical results presented in the case of the L-shaped membrane were obtained with a 40-digit 
precision in the eigenvalues, using the command N[,40] of Mathematical in this case, since I need to 
resort to larger grids I have worked with less digits precision using the command A^[] in Mathematica. 
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Figure 9. Energy of the ground state of the first isospectral membrane as a function 
of the number of grid points N. The horizontal hne is the precise result of [TU]. The 
set approaching the exact result from above (below) corresponds to e[^~^^ {e[^~^). 



The plot in Fig|9] displays the ground state energy of the first isospectral membrane 
calculated at different grid sizes. The horizontal line is the precise value given in flU\ . 
The set approaching this value from above (below) corresponds to the application of the 
method rejecting (accepting) the grid points falling on the border. The corresponding 
plot for the second isospectral membrane is almost identical and therefore it is not 
presented here. 

In Table |4] I report the energies of the first 30 states obtained using Richardson 
extrapolation [22] on the results for grids going from = 66 to = 120. The second 
and third columns are the energy of the first isospectral membranes obtained with the 
sets which reject (-Si^^) or accept {En ^) the grid points falling on the border, which 
as seen in the case of the L-shape membrane provide a sequence of numerical values 
approaching the exact eigenvalue from above and from below respectively. The last two 
columns report the analogous results for the second isospectral membrane. Notice that 
some of the energies in the third column are clearly incorrect. 

A further empirical verification of the isospectrality of the two membranes is 
presented in Figjlo| where I have plotted the asymmetry An = {En'^^ — En^^)/ {En'^'^ + 
E^^^) for the first 2000 states of the isospectral membranes. In this case 
is the energy of the n*^ state of the first (second) membrane obtained using Richardson 
extrapolation of the grids with = 114 and = 120. 



5. An unusual drum 



I will now consider a further example by looking at a particular membrane originally 
studied by Trott in [23]: this drum is shown in Fig, 12 and consists of a total of 
308 units squares which are joined into a rather irregular form. Theoretical and 
experimental studies carried out on drums with fractal or irregular boundaries have 
shown that the wave excitations for these drums are drastically altered [2^ |25| [26] : in 
particular, the Weyl law for these membranes is modified in a way which depends on 
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29 56976983 


29 56970152 


29 56975041 


29 56905778 


26 


31.48308074 


31.51241562 


31.48304984 


31.48393448 


27 


32.07624358 


32.16454642 


32.07622156 


32.08008665 


28 


32.21611001 


37.0118719 


32.21605287 


32.21393591 


29 


32.90535338 


27.9888228 


32.90537696 


32.90354978 


30 


34.13633502 


34.13929552 


34.13632946 


34.13632752 



Table 3. First 30 eigenvalues of the isospectral membranes obtained with Richardson 
extrapolation of the results obtained with grids from = 66 to iV = 120. 
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Figure 10. Left panel: logig of the asymmetry X = {E^n^'^ ~ E^n^'^) / {E^n^'^ + E^'n^^) 
for the first 2000 states of the isospectral membranes. En^^ {En^^ ) is the energy of 
the n*'' state of the first (second) membrane obtained using Richardson extrapolation 
of the grids with N = 114 and N ~ 120. Right panel: Blow-up of the previous plot 
for the first 100 states. 



the fractal dimension of the perimeter (see for example [2Z]), the so called Weyl-Berry- 
Lapidus conjecture. Recently the vibrations of a uniform membrane contained in a Koch 
snowflake have been studied in two papers, [281 EH] . 

The paper by Trott is both interesting in its physical and mathematical content and 
as an example of the excellent capabilities of Mathematica to handle heavy numerical 
calculations: as a matter of fact Trott uses a finite difference approximation of the 
Laplacian on a uniform grid and samples the membrane in 28521 internal points. 
Explicit numerical values for the first 24 modes are reported. 

I have therefore considered the same problem using the LSF with grids of different 
size (up to = 250 which leads to the same grid of [23j). Figure 13 displays the 
energy of the fundamental mode of this membrane as a function of the size of N . The 
dashed horizontal line in the plot represents the result of [23j, Ei = 6.64705: the points 
on the upper part of the plot correspond to N going from 50 to 250, with intervals of 
50. For these particular values of N the border of the membrane is sampled by the 
grid and therefore more accurate results are expected. The grid points on the border 
are rejected, which leads to eigenvalues which approach the exact results from above, 
as seen in the previous examples. The points in the lower part of the plot correspond 
to grid sizes varying from = 52 to = 148, excluding = 100: in this case the 
values approach the exact result from below, although in doing so they also oscillate 
reflecting the treatment of the border (a behaviour already observed in the case of the 
Africa membrane). As mentioned above the finest grid corresponds to sampling the 
membrane on 28521 internal points and therefore to working with a 28521 x 28521 
square matrix. Given that the matrix obtained with the LSF is a sparse symmetrix 
matrix, it is possible to deal efficiently with it in Mathematica, applying the Arnoldi 
method to extract a limited sequence of eigenvalues/eigenvectors. The reader will notice 
that in this example I have not considered the set corresponding to accepting the grid 
points falling on the border, as done in the case of the L-shaped and of the isospectral 
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Figure 11. Upper panel: Wave functions (absolute value) of the first isospectral 
membrane (ground state and 100*'' excited state); Lower panel: Wave functions 
(absolute value) of the second isospectral membrane (ground state and 100*'' excited 
state). A grid with iV 60 is used. 



membranes: although this set provides a sequence of values which uniformly approach 
the value at the continuum, the number of grid points sampled is quite large because of 
the large perimeter of the membrane. For example, for = 100, this set samples the 
membrane on 7029 points, compared with the = 3801 points used in the other set. 

The Figure also displays the improved ground state energies obtained using the 
"mesh refinement" procedure described in the Appendix (the three green points): the 
eigenvector for a given grid is extrapolated to a finer grid rejecting contributions in the 
"forbidden region" (i.e. falling outside the border of the membrane). The improved 
energy estimate corresponds to the expectation value of the Hamiltonian in this state 
and thus it requires no diagonalization. The results displayed in the figure correspond 
to extrapolation to a grid which is twice finer. 
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Figure 12. The unusual drum considered by Trott ^23^. The black area is the surface 
of the drum; the red points are the collocation points corresponding to iV = 50. 
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Figure 13. Energy of the ground state of the unusual drum as a function of N . The 
horizontal line is the result of |23j : the points approaching the horizontal line from 
above correspond to configurations where the border is sampled by the collocation 
points (and as discussed in the case of the L-shaped membrane are rejected). The 
green points correspond to the results obtained with the "mesh refinement" procedure 
described in the Appendix. 



6. Bound states in the continuum 



It is well known that the spectrum of the Laplacian with Dirichlet boundary conditions 
may contain bound states even for open geometries, in correspondence of crossings or 
bendings of the domain. For example, Schult et al.[3D] have studied the problem of 
two crossed wires, of infinite length, showing that such geometry supports exactly one 
bound state, localized at the crossing. Avishai and collaborators have also proved the 
existence of a bound state in the broken strip configuration for arbitrarily small angles, 
see (more recently Levin has proved the existence of one bound state in the broken 
strip for any angle of the strip [S2])- Goldstone and Jaffe have given a variational 




proof of the existence of a bound state for an infinite tube in two and three dimensions, 
provided that the tube is not straight. Other interesting configurations which support 
bound states in the continuum have been studied by Trefethen and Betcke 

The example which I will consider here is somehow related to the crossed wires 
configuration studied by Schult at al. I have considered a set of horizontal and vertical 
wires, of neglegible trasverse dimension, which are contained in a square box of size 
2. Calling n the number of wires in each dimension, is the number of crossings 
between these wires (for simplicity the wires are assumed to be equally spaced). This 
configuration can be easily studied in the present collocation approach, by sampling 
the wires on a grid and by then diagonalizing the Hamiltonian obtained following this 
procedure. The resulting energies calculated in this way will clearly depend on the 
spacing of the collocation grid, h, and diverge as h is sent to zero. To obtain finite 
results one needs to multiply these eigenvalues by /i^, which eliminates the divergence 
caused by the shrinking of the transverse dimension. Following this procedure I have 
studied different configurations, corresponding to choosing different value of n (going 
from 77, = 1 to n = 4) and I have found that a given configuration has precisely the same 
number of bound states as the number of crossing. These bound states happen to be 
almost exactly degenerate and correspond to wave functions which are localized on the 
vertices. 

In Table [4] I report the energy (multiplied by h"^) of the bound states and of the 
first unbound state (Egap) for the different configurations. These results have been 
obtained using a fine grid corresponding to h = 1/300 and show that the bound states 
are precisely as anticipated and they are essentially degenerate; the energy of the 
bound states and of the gap are also found to be almost insensitive to n, which can be 
interpreted as a sign of confinement of a state to the crossings. I have also checked the 
dependence of these results upon N (or equivalently upon h) observing that the energies 
can be fitted excellently as E = a + b/N'^; for example in the case of the ground state 
of the configuration with n = 4 I have obtained: E = 2.59874 — 44.6364/iV^. 
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Table 4. Energies of the bound states for configurations with different number of 
crossings, using N = 600, corresponding to a spacing h = 1/300. 



In Fig, 14 I have plotted the wave function of the ground state of the configuration 
corresponding to fi = 4 using a grid with = 500. The wave function is clearly localized 
at the crossings between the wires. Similar behaviour is observed for the remaining 15 
bound states. 



7. Collocation with conformal mapping 



The examples considered in the previous Sections show that it is possible to obtain the 
spectrum of the negative Laplacian over regions of arbitrary shape by using a collocation 
scheme, where the boundary conditions need not to be explicitly enforced on the border. 
Clearly, the precision of this approach should improve if the boundary conditions would 
be enforced exactly on the border of the membrane. One way of achieving this result is 
by mapping conformally the shape of the membrane into a square (or a rectangle), on 
whose border the LSF obey Dirichlet boundary conditions. I will discuss explicitly two 
examples of how this is done. 
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7.1. Circular membrane 

As a first example I consider a circular homogeneous membrane, which is exactly solvable 
(see for example [I]) and therefore it can be a useful tool to test the precision of the 
present method. 
The function 

/(z) = e"^ sn (^z F (^sin-i (^e"Tj _ i j _ i j (17) 

maps the unit square in the w complex plane into the unit circle in the complex z plane, 
as seen in Fig. [TSj Under this mapping the original equation, 

- AV^(w) = \i){w) (18) 
with Dirichlet boundary conditions on the unit circle, is mapped to 

- Ax{z) = \a{z)x{z) (19) 
with Dirichlet boundary conditions on the unit square. Here a{z) = |^|^ and 



eq. (19) describes the vibrations of a non-uniform square membrane. Although in the 
previous Sections I have restricted the application of the method to the case of uniform 
membranes of arbitrary shapes, the method can be applied also to inhomogenous 
membranes straightforwardly. Let me briefly mention how this is done. As a first 



step eq. ( 19 ) may be written in the equivalent form 
1 



criz 



Ax{z) = Xx{z) . (20) 



The operator O = -frA is evaluated on a uniform grid in the 2;-plane using the 
Little Sine Functions (LSF). The action of the operator over a product of sine functions 
can be calculated very easily, as explained in the previous Sections. To make the 
discussion simpler, I restrict to the equivalent one dimensional operator and make it 
act over a single LSF: 

Id? 1 

;Sk{h,N,x) = - ^—f — ^c'■^l'sj{h,N,x)sl{h,N,x) 



a{x) dx"^ ' ' ^-^ a{xj^ 

\- 1 (2) 



cl'>s,ih,N,x) . (21) 



The matrix representation of the operator over the grid may now be read explicitly 
from the expression above. The reader should notice that the matrix will not be 
symmetric unless the density is constant [|[| 

Using this approach I have considered grids with = 10, 20, . . . , 80 and I have 
have calculated the first four even-even eigenvalues, which are shown in Table |5j Taking 
into account the symmetry of problem I have used symmetrized LSF, which obey mixed 
boundary conditions (Dirichlet at one end and Neumann at the other hand): in this 
way, for a given value of the A^ a grid of (A^/2)^ points is used. As mentioned before the 

II In general the calculation of the eigenvalues and eigenvectors of non-symmetric matrices is 
computationally more demanding than for symmetric matrices of equal dimension. 




Figure 15. Unit square in the z plane and the corresponding unit cirle in the w plane 



reached through the trasformation (17) 



N 




El 


E2 


E3 




^4 


10 


5.7^ 


^5633618 


26.46056162 


30.55061880 


57.8f 


S187288 


20 


5.7^ 


^3347847 


26.37986506 


30.47598468 


57.60026669 


30 


5.7^ 


^3218252 


26.37564237 


30.47217988 


57.5J 


B626207 


40 


5.7^ 


^3196213 


26.37493961 


30.47155075 


57.5; 


B397911 


50 


5.7^ 


^3190167 


26.37474851 


30.47138009 


57.5f 


S336363 


60 


5.7^ 


^3187992 


26.37468004 


30.47131902 


57.5f 


S314408 


70 


5.7^ 


^3187059 


26.37465074 


30.47129291 


57.5f 


S305035 


80 


5.7^ 


^3186606 


26.37463653 


30.47128023 


57.5J 


B300497 


LSQs 


5.7^ 


^3185971 


26.37461646 


30.47126209 


57.5J 


B294087 


Exact 


5.7^ 


^3185962 


26.37461642 


30.47126234 


57.5J 


B294090 



Table 5. Even-even spectrum of the circular membrane: first four eigenvalues 



exact eigenvalues for this problem are known (the zeroes of the Bessel functions): these 
are reported in the last row. 

In fig. 17 I have plotted the lowest eigenvalue of the circular membrane 
corresponding to different and I have fitted these points using functions like cq+ci/N'^, 
with r = 3,4,5 (the dashed, solid and dotted lines in the plot). This plot shows that 
the leading (non-constant) behaviour of the numerical energy for ^ 1 is 
Taking into account this behaviour I have considered the quantity 

2 

Or 



8 

E 

k=l 



«1 



(lOA;)"^ 



(22) 



where Q = 8 and I have obtained the coefficients by minimizing Hq (notice that 
this expression takes into account the leading behaviour just discussed). The row 
marked as LSQs displays the quite precise results obtained following this procedure. 
I would like to discuss briefly a different issue. In [M] Gottlieb has used the Moebius 
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Figure 16. Density of the inhomogeneous square membrane isospectral to the 
homogeneous circular membrane. 




Figure 17. Energy of the ground state of the circular membrane. The dashed, solid 
and dotted lines correspond to fits using functions like cq + ci/N^', with r = 3,4,5 
respectively. 



transformation 

fg{z) = {z- a)/(l - az) (23) 

to map the unit circle onto itself. This mapping transforms the homogenoeous Helmoltz 
equation for a circular membrane into the inhomogeneous Helmoltz equation for a 
circular membrane with density 

Pix, y) = If'M" = ^^TTT^^^tW • (2^^ 

[(1 — ax)'^ + a^y^l 

Gottlieb uses this result to conclude that membranes corresponding to different 
densities, i.e. different values of a, are isospectral, thus providing a negative answer to 
the famous question "Can one hear the shape of a drum?" , posed by Kac in [TB] . I wish 
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to move our discussion on computational grounds: for a given a the mapping of eq. (23) 



deforms the grid inside the unit circle; as a is changed, the grid points move, as shown 



in Fig. 18 The case a = is plotted in the right panel of Fig. 15 Clearly, if the density 



of the membrane is constant, or symmetric with respect to the center, one expects that 



a = provide the best grid. In Fig. 19 I have plotted the logarithm of the difference 
between the approximate and exact energy for the ground state of a circular membrane, 
A = LogiQ{E]\r — E^xact), using three values of a (a = 0, 0.4 and 0.8). These numerical 
results confirm the prediction made: stated in different terms one can conclude that for 
a given problem one can improve the numerical accuracy of a calculation by selecting an 
optimal grid among those obtained through a conformal map of the region onto itself. 
The optimization of the parameter a depending on the specific problem considered is 
in the same spirit of the variational approach used in [TH| [T3] and could provide a 
useful computational tool to boost the precision of the results. 





Figure 18. Grid obtained witli tlie Moebius map corresponding to a = 0.5 (left) and 
a = -0.8 (right). 




Figure 19. A = LogiQ{EN — E^xact) using three values of a (a = 0, 0.4 and 0.8 from 
bottom to top). 
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7.2. Circular waveguide 

The second example of application of conformal mapping to the solution of the 
Helmholtz equation is taken from the paper of Kuttler and Sigillito p] (this problem 
was also studied earlier by Moler, in ref. [101] of [2]). 



In Figj20]two regions of the plane are displayed: the left plot corresponds to a square 
of side TT centered on the origin in the z = x + iy plane; the right plot corresponds to a 
circular waveguide with circular ridges in the w = u + iv plane. The function w = tan | 
maps the first region into the second one. 

As I have shown for the case of the circular membrane, the homogeneous Helmoholtz 
equation over the second region may be transformed into an inhomogeneous Helmholtz 
equation over the square: 

- AU{z) = Xa{z)U{z) . (25) 

In the present case a{z) = |^|^ = (cosx + coshy)^ and Dirichlet boundary conditions 
are assumed on the borders of the two regions. 

In Tables 1,2 and 3 of their paper, Kuttler and Sigillito report different estimates for 
the first 12 even-even eigenvalues, obtained using different approaches. In Table 2 they 
also apply Richardson extrapolation to the eigenvalues obtained with finite difference. 
In the case of the ground state of this membrane they also mention the precise value 
obtained by Moler using the method of point matching 

Ai = 7.5695769 (26) 



In Table 6 I report the even-even eigenvalues of the eq. (25) obtained using 



collocation with different values of A^. The results corresponding to the ground state are 



plotted in Fig. 21 and fitted using functions like cq + ci/N^, with r = 3, 4, 5 (the dashed, 
solid and dotted lines in the plot). This plot proves that the leading (non-constant) 
behaviour of the numerical energy for ^ 1 is l/N'^, as for the circular membrane. 
The results in the Table have also been extrapolated using a least square approach 



2 

Ear, ' 



(lOA;)"- 



(27) 



k=l 

where Q = 7,8 and a„ are coefficients which are obtained by minimizing Sg. Notice that 
this expression takes into account the leading l/N"^ behaviour just discussed. The rows 
marked as LSQy^s display the results obtained following this procedure (the comparison 
between the results for Q = 7 and Q = 8 gives an indication over the precision reached): 
in particular the energy of the ground state reproduces all the digits of the result 
obtained by Moler. It is also remarkable that the energies obtained with the conformal- 
collocation method decrease monotonically when the number of collocation points is 
increased (the only exception is represented by the £"10 for A^ = 10, probably due to the 
limited number of collocation points). 

As a technical remark, one should notice that the results corresponding to a given 
value of A^ are obtained using a set of A^/2 symmetric (even) functions for each direction. 
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Table 6. Even-even eigenvalues of the problem of eq. (25 1 using collocation with Little 
Sine Functions (LSF). 
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thus reducing the computation load by a factor of 4. The results displayed in this table 
should be compared with the analogous results of Table 2 of [2], which were obtained 
using finite difference. 



y 




Figure 21. Energy of the ground state of the circular waveguide. The dashed, sohd 
and dotted hues correspond to fits using functions hke cg + ci/iV, with r ~ 3,4,5 
respectively. 

8. Conclusions 

In this paper I have used a collocation method based on LSF to obtain the numerical 
solutions of the Helmholtz equation over two-dimensional regions of arbitrary shape. A 
large number of examples has been studied, illustrating the great potentialities of the 
present method. Among the principal virtues of this method I would like to mention its 
generality (it can be applied to membranes of arbitrary shapes, including inhomogeneous 
membranes, and to the Schrodinger equation - although I have not done this in the 
present paper), its simplicity (the matrix representation of the Helmholtz operator is 
obtained directly by collocation, and therefore it does not require the calculation of 
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Figure 22. Upper panel: Even-even wave functions (absolute value): ground state and 
100*'' excited state of the circular waveguide; Lower panel: Even-even wave functions 
(absolute value): 200*'* and 300*'* excited states of the circular waveguide. A grid with 
iV = 80 is used. 

integrals) and the possibility of combining it with a conformal mapping, as done in 
the last Section. In this last case, a rapid convergence to the exact eigenvalues is 
observed as the number of grid points is increased. In the case where the border is 
not treated exactly it has also been observed that the method provides monotonous 
sequences of approximations to the exact eigenvalue either from above or from below. 
Readers interested to looking at more examples of application of this method may find 
useful to check the gallery of images which can be found at 

\protect\vruleuwidthOpt\protect\href {http : //f e j er . ucol . mx/paolo/ druin}{http : //f e j er . ucol . 
Appendix A. Mesh refinement 

Although the collocation method described in this paper allows one to obtain precise 
solutions to the Helmholtz equation over domains of arbitrary shape, in general the 
Dirichlet boundary conditions are not enforced exactly over all the boundary. As 
discussed in Section [7] the best approach consists of introducing a conformal map, which 
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allows one to go from the original problem to an inhomogenous Helmoltz problem over 
a square: in such case the Dirichlet boundary conditions are imposed exactly and rapid 
convergence to the exact solutions is observed. In general, however, finding such a 
conformal map can be a difficult task and therefore the first approach may be more 
appealing. I will discuss here a simple procedure to "refine" the results obtained by 
direct collocation of the Helmholtz equation over the grid. The fundamental observation 
is that the LSF that we have used do vanish on the grid points on the border and external 
to the membrane, but they are nonzero in all the other points external to the membrane. 
Therefore the cumulative effects of the LSF internal to the membrane can be seen also 
outside the membrane, although it will tend to disappear as the number of grid points 
is increased. This solution, to increase the number of grid points, may be the most 
obvious but it is certainly not appealing computationally, since increasing the number 
of grid points strongly increases the computational cost (remember that the number of 
matrix elements grows as A^"^). However we can use much simpler procedure, which does 
not require any additional diagonalization. Call N the parameter defining the size of 
the grid: a point in this grid is described by the direct product of the LSF in the x and 
y directions. In the Dirac notation we write 

{x,y\k,k')h Sk{h,x)sk'{h,y) , (A.l) 

assuming for simplicity that the grid has the same spacing in both directions. Let us 
now concentrate on one of the LSF, say the one in the x direction: we take a finer grid, 
with a spacing h' — h/l, where I is a integer. The new grid contains now {IN — 1) points, 
including obviously the original grid points. However, it is clear that the original LSF 
can be decomposed in the new grid as 

IN /2-l 

Sk{h,x)^ ^ Sk{h,Xj) Sj{h/l,x), (A.2) 

j=-lN/2+l 

where Xj = 2Lj/{lN) are the new grid points. Notice that this relation is exact. 

The wave function of the n'^^ state obtained from the diagonalization of the 
[N — 1) X [N — 1) hamiltonian reads 

ipn{x,y) = ^'vp Sk{K){h,x) Sk'{K){h,y) 

^ K 

^ IN /2-1 /7V/2-1 

^h^^K^ ^ Sk{K){h,Xj) Sj{h/l,x) Sk'(K){h,yf) Sf{h/l,y) , 

where the n*'* eigenvector. Clearly '4>n{x,y) differs from even in points of the 

refined grid which fall outside the membrane profile. We introduce a new matrix whose 
elements are given by 

_ J i/ ^ B 

^-■'-1 1 i/(x,,y,0 e B ^""-'^ 
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and rewrite the wave function "purged" on the refined grid as 

IN /2~1 ZAf/2-1 



AT 



j=-lN/2+l j'=-lN/2+l 
in) 



where 

K 

and A/" is a normahzation constant that ensures that 
i)l{x,y)dx dy = l. 



It is easy to show that 



To simphiy the notation I define: 
AT.-. 



^33' = y^/ 



and thus write: 



/Af/2-1 /Af/2-1 



(A.4) 



(A.5) 



(A.6) 



(A.7) 



j=-lN/2+l j'=-lN/2+l 

On the other hand we may also calculate the expectation value of the Hamiltonian 
in this state 



{H)n = - / i^n{x,y) AiJri{x,y) dx dy 



- E 



Vjj'Vrr' 



Sj{h/l,x) Sf{h/l,y)ss{h/l,x) Ss'{h/l,y) 



jj'rr' 



(A.8) 



33 r 



where C(2) is the matrix for the second derivative on the refined grid. An example of 



application of this procedure is shown in Fig, 13 
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